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Abstract 

We have performed quantum molecular-dynamics simulations for methane under shock compres- 
sions up to 80 GPa. We obtain good agreement with available experimental data for the principal 
Hugoniot, derived from the equation of state. A systematic study of the optical conductivity spec- 
tra, one-particle density of states, and the distributions of the electronic charge over supercell at 
Hugoniot points shows that the transition of shocked methane to a metallic state takes place close 
to the density at which methane dissociates significantly into molecular hydrogen and some long 
alkane chains. We predict the chemical picture of the shocked methane with respect to the pair 
correlation function . In contrast to usual assumptions used for high pressure modeling of methane, 
we find that no diamond-like configurations occurs for the whole density-temperature range studied. 
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I. INTRODUCTION 



The nature of methane under extreme conditions has recently drawn extensive attentions 
due to the ongoing scientific and technological interest. As methane is one major constituent 
of the "ice" layers in Uranus and Neptune, where the pressure ranges between 20 and 600 GPa 
and temperature between 2000 to 8000 K |l, 2|, its properties are of considerable impact on 
the internal evolution and energetics of these giant planets. The knowledge of its equation 
of state and electrical conductivity at these conditions are necessary both for theoretical 
modeling of the composition of these planets and for understanding its contribution to the 
planetary magnetic fields, which are caused by convective dynamo action of electrically 
conducting fiuid jsl- Thus, it is desirable to determine the transformations in methane 
induced by nonequilibrium phenomena such as shock waves and detonations. 

A lot of efforts have been devoted to exploring the methane chemical dissociation under 
extreme conditions. At high static pressures between 10 and 50 GPa and temperatures of 



about 1000 to 3000 K, diamond-anvi 
down to form diamond and hydrogen 



cell experiment has suggested that methane breaks 
4 1, in qualitative agreement with Ree's prediction j^. 



However, on the methane Hugoniot [6J different results were reported. No dissociation was 
found in shocked methane for compressions up to 42 GPa [7], while Nellis et al. observed the 
increase of electrical conductivity in methane shocked by two-stage light-gas gun, which was 
attributed to the decomposition into a substantial concentration of molecular hydrogen at 
36 GPa and 3900 K [8]. The discrepancy may be caused by the shorter time scale of shock 
wave experiments. Concerning the theoretical side, empirical potentials jo], [lo| and tight- 
binding simulations |ll| have been used to conduct molecular dynamics simulations of the 
shock compression of methane, whi 
molecular dynamics simulations 
forming of chemical bonds. FPMD simulations with 16 methane molecules for 2 picoseconds 
show that some polymerization reactions occur at 100 GPa and 4000 K and that diamond 



e the use of density functional theory in first-principles 
isl provides more accurate modeling of the breaking and 



12| . Recent ab initio evolutionary simulations 



brmation take place only above 300 GPa 

14j | have systematically investigated the phase diagram of methane under pressure, and 
confirmed the dissociation of methane at high pressures. Whereas, the anharmonism was 
neglected in calculations, which may affects the stability of methane. 

To date, although a number of explanatory and predictive results in some cases have al- 
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ready been provided by experimental and theoretical studies, many fundamental questions 
of methane at extreme conditions are still much controversial and under debate. The system- 
atic study of equation of state, especially the Hugoniot curve under extreme conditions is far 
lacking now, which is essential in this context. Furthermore, in adiabatic compressions, the 
nonmetal-metal transition of hydrocarbons has been a major issue recently. For example, 
benzene jl5| has been reported to transform metallic phase when decomposing into hydro- 
gen under extreme conditions. To date there are still no such data explored for methane, 
especially due to the difficulties in applying current opacity models at these particular condi- 
tions. It is thus highly needed to study the electronic structure, Hugoniot EOS, and optical 
properties of shocked methane using an efficient method consistently. Quantum molecular 
dynamics (QMD), where the active electrons are treated in a full quantum mechanical way 
within the finite temperature density functional theory (FT-DFT), has been proven a suc- 



cess: 
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ul tool to calculate the properties of complex plasmas under such extreme conditions 



2l| . When combined with the Kubo-Greenwood formulation, the method produces a 



consistent set of material, electrical, and optical properties from the same simulation and 
can be applied to various systems. 

In the present work, we perform QMD simulations for the system with densities and 
temperatures ranging from 0.6 g/cm^ and 175 K to 1.5 g/cm^ and 7000 K along the prin- 
cipal Hugoniot on an intermediate time scale. The analysis of the concentration of molec- 
ular species along the Hugoniot based on pair correlation function (PCF) indicates that 
methane decomposes into molecular hydrogen and saturated hydrocarbons first, and then 
into long alkane chains at higher pressures and temperatures, which consequently leads 
to the nonmetal-metal transition. The dynamic conductivity, calculated using the Kubo- 
Greenwood formula, along with electron density of states and charge distribution, confirms 
the occurrence of nonmetal-metal transition around 55 GPa. This paper is organized as fol- 
lows. The simulation details are briefiy described in Sec. II; the PCF, which is used to study 
the dissociation of methane, and Hugoniot curve are given in Sec. Ill; in Sec IV nonmetal- 
metal transition, optical properties and electronic properties are discussed. Finally, we close 
our paper with a summary of our main results. 
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II. COMPUTATIONAL METHOD 



In the present application, the molecular dynamics trajectories are calculated by employ- 
ing the Vienna ab initio simulation package (VASP) plane-wave pseudopotential code devel- 



oped at the Technical University of Vienna 
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23|. In these simulations, the electrons receive 



a fully quantum mechanical treatment through solving the Kohn-Sham (KS ) eq uations for 



a set of orbitals and energies within a plane-wave, FT-DFT formulation 



25|, where the 



electronic states are populated according to the Fermi-Dirac distribution at temperature 



Tg. The all-electron projector augmented wave (PAW) method 26|, |27| is adopted and the 



exchange-correlation energy is described employing the Perdew-Wang 91 parametrization of 
the generalized gradient approximation (GGA) 28|. Atoms move classically according to 
the forces, which originate from the interactions of ions and electrons. 

We perform finite-temperature, fixed-volume molecular dynamics simulations for selected 
densities ranging from 0.6 to 1.5 g/cm^ and temperatures from 175 to 7000 K that highlight 
the single-shock Hugoniot region. We use 27 carbon and 108 hydrogen atoms (twenty-seven 
methane molecules) in a cubic cell and fix the plane-wave cutoff at 600.0 eV which is tested 
to give good convergence for both total energy and pressure. The Brillouin zone is sampled 
with the F point for molecular dynamics and 4x4x4 Monkhorst-Pack 29| scheme k points 



for the electronic structure calculations. Integration of the equations of motion proceeds with 
time steps of 0.5-1.0 fs for different pressure-temperature ranges. Typical simulations run 
for 8000-16000 steps with the time scale about 8 ps. We let the system equilibrate for 4000- 
8000 steps and calculate properties using the final 4000-8000 steps. The isokinetic ensemble 
(NVT) is employed for the ions, where the ion temperature Tj is fixed using velocity scaling. 
The electron temperature Te is in turn set to that of the ions Tj based on the assumption 
of local thermodynamical equilibrium. The accuracy of our calculations is examined by the 
bond length of CH4 molecule in its ground state and the result is 1.09, which agrees well 
with the experiment. 

III. SHOCK EOS AND PCF 

A precise description of material properties, such as EOS, is important to the accurate 
calculation of electrical and optical properties. A crucial measure for theoretical EOS data 
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is the principal Hugoniot curve. It describes the shock adiabat through a relation between 
the initial and final internal energy, pressure and volume, respectively, {Eq, Pq, Vq) and 
{El, Pi, Vi) according to the following Rankine- Hugoniot equation jsol : 



{El -Eo) + ^ {Vi - Vo) (Po + Pi) = 0, (1) 

where the internal energy E equals to the sum of the ion kinetic energy |A;s7i, the time 
average of the DFT potential energy and zero-point energy with ks the Boltzmann constant. 
The pressure consists of contributions from the electronic Pe and ionic Pj components, which 
come from, respectively, the derivatives taken with respect to the KS electronic orbitals and 
the ideal gas expression since ions move classically. We thus have P = Pe + pn/c^T, where 
Pn is the number density. For the methane principal Hugoniot curve, the initial density is 
Po = 0.423g/cm^ and the initial internal energy Eq = 23.30 eV/molecule at a temperature of 
111 K. The initial pressure can be neglected compared to the high pressure of the final state. 
The Hugoniot points are determined in the following way. For a given density, the periodic 
crystalline sample of the corresponding size was first constructed. A fixed volume cubic 
supercell of 27 carbon and 108 hydrogen atoms (27 methane molecules), which is repeated 
periodically throughout the space, forms the elements of the calculation. A series of quantum 
molecular dynamics simulations are performed for different temperatures T. Following it, 
the internal energies and pressures are determined correspondingly and then fitted to a cubic 
function of T. For a given density and a set of temperatures, we plotted {Vq — Vi) (Pq + Pi)/2 
along with {Ei — Eq) as a function of temperature. The intersection of the two curves fixes 
the principal Hugoniot point {Ei, Pi, Vi) that satisfies Eq.(l). The particle velocity Up and 



shock velocity Ug are then determined from the other two Rankine- Hugoniot equations |30| . 

= Vq [1 — {up/ug)] and Pi — Pq = pQUgUp. The principal Hugoniot points of methane 
derived from Eq. (1) are listed in Table I. 

We display our simulated Hugoniot curve for methane in Fig. 1, along with experimental 
data and results from previous theoretical works. From Fig. 1 we find very good agreement 
with the experimental results all along the single-shock Hugoniot, while the data from the 
empirical potential simulations based on the adaptive intermolecular reactive empirical bond 
order force field presents prominent discrepancy with all the other results, which may be 
caused by the inaccurate modeling of the breaking and forming of chemical bonds. Our 
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FIG. 1: Principal Hugoniot of methane. For comparison, the present QMD results, together with 
the previous experimental data (Ref. 6 and Ref. 7) and other theoretical results (Ref. 10 and Ref. 
13) are all provided. 

simulated pressures and temperatures using QMD show a systematic behavior along the 
Hugoniot curve, except for the region with a small break in slope of temperature plot between 
1.16 and 1.26 g/cm^, which usually implies the dissociation of molecular species in the media 
and thus suggests that methane is dissociating in this density region. 

To clarify the structural change in methane under shock conditions, we calculate the 
PCFs for each pair of atom types, which give the possibility of finding an atom of a given 
type at a given distance from a reference atom. The PCFs, together with atomic structure 
along the principal Hugoniot of methane, are presented in Fig. 2. At the lowest density of 
p = 0.423 g/cm^, the PCF of C-H bond qc-r (r) peaks at about I.OQA, which corresponds to 
the equilibrium internuclear distance of the C-H bond in methane molecule. At this density, 
the hydrogen correlation function does not show a maximum at a distance corresponding to 
the equilibrium distance of hydrogen molecule Th-h = 0.75 A. Meanwhile, no peaks occur in 
the PCF of C-C bond at 1.50A or 1.54A, which corresponds to the typical C-C bond length 
in diamond or saturated hydrocarbons. Therefore, the PCFs at 0.423 g/cm^ indicate that 
methane remains its ideal molecular configurations without dissociation. As the density is 
compressed to p = 1.16 g/cm^, methane molecules dissociate with the increase of density 
and temperature, which is indicated in Fig. 2(b) by the significant reduction and broadening 
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TABLE I: Points along the principal methane Hugoniot derived from DFT-MD simulations at a 
series of density (p), pressure (P), and temperature (T). 



p 


P 


T 


(g/cm^) 


(GPa) 


(K) 


0.60 


2.06 


176 


0.73 


5.02 


330 


0.86 


11.89 


1441 


1.05 


25.42 


2614 


1.16 


36.74 


3420 


1.26 


45.55 


3620 


1.40 


60.20 


4509 


1.50 


84.00 


6875 



of the maxima of g'c-H {t) around I.OQA. On the other hand, the g^-n and gc-c PCFs at this 
density indicate that small amount of hydrogen and saturated hydrocarbons form molecules 
upon dissociation of methane. Consistently, the atomic configuration consists of a mixture 
of methane and ethane with the residual hydrogen atoms both in molecular and atomic 
forms as shown in the inset of Fig. 2(b). As the density is further increased to 1.26 g/cm^, 
maximum of the PCF gc-n {f) continues to reduce and broaden while the peaks of g^n-H {f) 
and gc-c {f) increase, which indicates that methane further dissociates into hydrogen and 
some hydrocarbons of higher molecular weight. This trend persists up to a density of 1.4 
g/cm^, except that a new species of long alkane chains show up in the atomic configuration. 
It is indicative of the formation of C-C bonds favored with increasing density, though still no 
diamond-like carbons form in the density range considered here, in agreement with earlier 
predictions. 

IV. DYNAMIC OPTICAL AND ELECTRONIC PROPERTIES 

Following the QMD simulations, a total number of 10 to 15 configurations arc selected 
from an equilibrated (in an average sense) portion of the molecular dynamics run, typically 
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FIG. 2: Pair-correlation functions for C-C (red or gray line), C-H (black line), H-H (green or 
light gray line) along the principal methane Hugoniot. The atomic structure, where carbon and 
hydrogen atoms are denoted by gray and white balls, respectively, is also provided in the insets. 
(a)p = 0.423 g/cm3, T = lllK; {h)p = LlGg/cm^, T = 3420K; {c)p = 1.26g/cm3, T = 3620K; 
(d)/9 = 1.40g/cm3, T = 4509 K. 

sampling the final picosecond of evolution. The configurations are spaced at time steps sep- 
arated by at least the correlation time, the e— folding time of the velocity autocorrelation 
function. The calculated properties should be averaged over the number of representative 
configurations or snapshots employed. For each of these configurations, the electrical con- 
ductivity is calculated using the most general formulation given by the Kubo-Greenwood 
formulation, without particular assumptions made on the ionic structure or on the electron- 
ion interactions. In the framework of the quasi-independent particle approximation, the 
Kubo-Greenwood formulation 
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32] gives the real part of the electrical conductivity as a 



function of frequency cj. 



„ N N 1, 

= ^E^(k)EEE[/(^-k)-/(e„k)] 

k j=l 1=1 a=l 

|2 



X |(^j,k| V«|*i,k)r5(ej,k-ei,k-^), (2) 
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where uj is the frequency, is the volume of the supercell, N is the total number of energy 
bands used, ^i,k and e^^k are the electronic eigenstates and eigenvalues for the electronic state 
i at k, / (ej, k) stands for the Fermi distribution function, and w (k) represents the k-point 
weighting factor. Other properties can be directly derived from the frequency-dependent 
real part of the electrical conductivity. An application of the Kramer-Kronig relation yields 
the imaginary part 02 (<^) as 

-^M = -;^/(^<'-. (3) 

where P stands for the principal value of the integral. The real and imaginary parts of 
dielectric function, in turn, follow immediately from the electrical conductivity as 

eiH = 1 (T2(a;), (4) 



47r 

62 (u;) = — (7i (a;) . (5) 

UJ 

And then the real n [fJ] and imaginary k (a;) parts of the optical refraction index have 
relation with the dielectric function by a simple formula. 



e {uj) = ei {uj) + ie^ {uj) = \n {uj) + ik (u;)]^ . (6) 

Finally, the reflectivity r {uj) and absorption coefficient a {uj) can be determined from these 
quantities. 



[1 + n (cujj + k {uj) 



4:71 

a (c) = (c) . (8) 

In Fig. 3 the behavior of the frequency-dependent conductivity ai {uj) at different Hugo- 
niot points is reported in "raw" data form. It is found that ai {uj) at these four different 
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FIG. 3: The real part of the dynamic electrical conductivity along the principal Hugoniot. Data 
have been averaged over 10 to 15 uncorrelated MD configurations. 

densities exhibit a uniform feature that they peak around 10.0 eV, which can be associated 
with transitions to the lowest excited states, and vanish above 25.0 eV. With the increase 
of density and temperature along the principal Hugoniot, the shape of cti (uj) mostly keeps 
the same, but the main peak moves to lower frequencies, which consequently leads to a 
significant increase in dc conductivity, given as adc = Hmai (u). This significant change of 
the dc conductivity along principal Hugoniot is displayed in Fig. 4. For comparison, the 
experimental data is also included. The difference between our calculated results and the 
experimental results may comes from the under-prediction of the Hugoniot temperatures 
for methane by quantum molecular dynamics. The dc conductivity is negligible for pres- 
sures up to 40 GPa, where methane exhibits insulating behavior. Then, it increases to 128 
fi~^cm~^ for pressures between 40 and 50 GPa. For pressures above 55 GPa, the dc conduc- 
tivity increases to a value lager than 1000 Q~^cm^^. A decisive physical assertion to define 



the metal 
by Mott 



icity of such disordered system has been proposed by loffe et al. and developed 
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34j |. which states that any high-temperature, disordered system will remain 
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FIG. 4: The dc conductivity of our results and data from Ref. 8 are plotted along the principal 
Hugoniot. 

metallic or indeed attain metallic status if the characteristic mean free path of the valence 
electrons exceeds the mean distance between the constituent atoms or molecules providing 
those carriers of electrical current. This simple but powerful argument leads to an estimate 
of the minimum metallic conductivity of about 2000 Q^^cm~^ for fluid hydrogen, rubidium, 
caesium and mercury, whereby transitions from non-metallic to metallic fluids have recently 
been observed in the nominally non-mitallic chemical elements hydrogen, oxygen and ni- 



trogen at high pressures and temperatures 
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With the similar minimum metallic 



conductivity for fluid methane, it is indicated that nonmetal-metal transition takes place in 
shocked methane. This nonmetal-metal transition clearly suggests that the dissociation of 
the molecular fluid along the Hugoniot has major consequences to the electrical conductivity 
of the system. Interest ing ly, such nonmetal-metal transition in fluid hydrogen occurs under 



similar conditions 
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Pressure (GPa) 

FIG. 5: Optical reflectivity of shocked methane for wavelengths of 404 (black square), 808 (red 
circle) and 1064 nm (blue triangle) along the Hugoniot. 

As the optical reflectivity of fluid has been successfully probed in dynamic compression 
experiments (39|, we show in Fig. 5 the reflectivity r (lu) at typical wavelengths of 404, 808 
and 1064 nm spanning the visible spectrum. It can be seen that a measurable reflectivity 
arises from 0.05 to 0.35-0.68 for the principal shocks, which is related with the high-pressure 
nonmetal-metal transition. For pressures between 25 and 45 GPa reflectivity increases mono- 
tonically and is essentially equal for all three frequencies, while around 55 GPa the increase 
becomes shaper due to the nonmetal-metal transition, after which the increase slows down. 

Furthermore, we examine the variation of the electron density of states (EDOS) along the 
Hugoniot to clarify the origin of this nonmetal-metal transition. Still use the selected ten 
to flfteen conflgurations used in the above optical properties calculations, we calculate the 
electron density of states. Each of these conflgurations repeats periodically throughout the 
space and forms the elements of the calculation. In Fig. 6 the averaged EDOS is presented 
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FIG. 6: Electronic density of states of liquid methane along the principal Hugoniot: {a)p = 
1.05g/cm^ T = 2614K; {h)p = l.lGg/cm^, T = 3420K; (c)p = 1.26g/cm^ T = 3620K; 
(d)p = 1.40g/cm3, T = 4509K. Data have been averaged over 10 to 15 uncorrelated MD con- 
figurations. The zero of the energy scale shows the position of the Fermi level. 

for different Hugoniot points. As can be seen, at p = 1.16 g/cm^ and T = 3420 K, EDOS 
clearly exhibits a gap around 2.4 eV, where only thermally activated electron transport 
occurs as in semiconductors. With the density and temperature increasing, the band gap is 
reduced gradually and closed eventually. And thus a higher, metal-like conductivity follows. 

Another way to characterize this behavior is the change of the charge density with in- 
creasing density as shown in Fig. 7, from which it is clearly seen that a transition from 
a disconnected network (left panel) to a connected one (right) occurs as density increases. 
This suggests that the electrons are still localized for the lower density, while transient fila- 
ments and lusters form as percolating path for the higher density, which leads to a metal-like 
behavior. 
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FIG. 7: Typical contour plots of the charge densities in units of e/A for 0.86 g/cm^ at 1441 K 
(left panel) and 1.4 g/cm^ at 4509 K (right panel) along a slice through the simulation box. 

V. CONCLUSIONS 

In summary, through systematic QMD simulations we have determined the shocked EOS 
and clarified the high-pressure nonmetal-metal transition of the fluid methane. The increase 
of electrical conductivity with pressure can be ascribed to the dissociation of the molecular 
fluid. However, no diamond-like carbon forms in the density range considered here. In 
addition, the optical properties of warm dense methane are also calculated, from which 
an experimentally measurable increase in the reflectivity associated with the high-pressure 
nonmetal-metal transition has been found. It is expected that our simulated results of 
shocked methane would have a strong influence on models for planetary interiors. 
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